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Abstract 

We apply inverse statistical-mechanical methods to find a simple family of optimized isotropic, 
monotonic pair potentials, under certain constraints, whose ground states for a wide range of 
pressures is the diamond crystal. These constraints include desirable phonon spectra and the 
widest possible pressure range for stability. We also ascertain the ground-state phase diagram for a 
specific optimized potential to show that other crystal structures arise for other pressures. Cooling 
disordered configurations interacting with our optimized potential to absolute zero frequently leads 
to the desired diamond crystal ground state, revealing that the capture basin for the global energy 
minimum is large and broad relative to the local energy minima basins. 
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Advances in the field of self-assembly, devising building blocks (e.g., nanoparticles and 
polymer chains) with specific interactions to form larger functioning materials, are proceed- 
ing rapidly and hold great promise to produce unique colloidal and polymer systems.-~- 
In the past several years, inverse statistical-mechanical methods have been formulated that 
yield optimized interactions that robustly and spontaneously lead to a targeted many-particle 
configuration with desirable or novel bulk properties.- This inverse approach provides a pow- 
erful and systematic means of directing self assembly with exquisite control. Recent studies 
have used inverse methods to find optimized isotropic (non-directional) interactions, subject 
to certain constraints, that yield novel targeted ground states, such as low-coordinated crys- 
tal structures.- This includes the three-fold coordinated honeycomb (or graphene) structure 
in two dimensions^ and the tetrahedrally-coordinated diamond crystal in three dimensions,- 
initial studies of which involved isotropic pair potentials with multiple wells. 

Are multiple wells required to achieve low-coordinated crystal ground states with isotropic 
pair interactions? We have recently shown that inverse statistical-mechanical techniques al- 
low one to produce robustly both the square lattice and honeycomb crystal in two dimensions 
via monotonic convex pair potentials.- 1 ^ In the present work, we employ inverse techniques 
to obtain a simple family of optimized isotropic, monotonic pair potentials whose ground 
states for a wide range of pressures is the diamond crystal. 

Using the forward approach,— it was established over a decade ago that the diamond 
crystal can be stabilized for a range of densities by an isotropic, monotonic pair potential 
devised to model star polymers.— These authors used free-energy calculations to find the 
phase diagram and validate their conclusions. Moreover, the potential possessed stable 
phonon spectra over the predicted ground-state parameter regime.— A forward approach 
was used in another study^ to examine only lattice energy sums at zero temperature for 
a relatively small set of Bravais and non-Bravais lattices for an isotropic, monotonic pair 
potential. It was found that the diamond crystal was stable for a certain pressure range. 
These authors recognized the limitations of this restricted investigation, which excluded 
both phonon spectra calculations and annealings to zero temperature from liquid-like initial 
conditions in order to validate that the diamond was indeed the ground state. 

Here we use a simpler functional form for a monotonic radial (isotropic) pair potential 
function v(r) that obeys certain important conditions on the second derivative with respect 
to the radial interparticle distance r established in Refs. ^ and [9). Specifically, we propose a 
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potential function of the form 

v{r)=e{l + a l r - + a 2 [^e-^\ (1) 

where e and a, respectively, define the energy and length units, and a\ and a 2 are dimen- 
sionless parameters. Equation ([1]) is chosen for its simplicity and because it allows for the 
desirable features of the second derivative described below. The potential function ([I]) is 
strictly convex for all r beyond a small cutoff distance for a large range of parameters. In 
this study, we restrict ourselves to such potentials that are convex for r > 0.1a.— 

We use an iterative two-step inverse procedure to determine the optimized parameters of 
the potential function ([1]) under certain constraints that yields the diamond ground state 
for a range of pressures. We define the dimensionless pressure p* and density p* as follows: 

p* = pa 3 /e, p* = a 3 p. (2) 

The first step of the optimization procedure involves choosing an initial set of "competitor" 
configurations. Then we determine the parameters a\ and a 2 that maximize the ratio be- 
tween the maximum and minimum pressures p ma x and p m in for which the diamond crystal 
has a lower enthalpy than any competitor configuration, since we are working in the isobaric 
ensemble. The second step involves a rapid cooling procedure within a simulation box under 
periodic conditions in the isobaric ensemble, implying that the simulation box is deforming 
and changing volume. We start by choosing initial lattice vectors that define the simulation 
box within which there is a Poisson distributed set of N particles interacting with the po- 
tential function (CQ) and the parameters from the first step. Cooling is then carried out using 
a quasi-Newton method, which has similar basins of attraction as obtained from steepest- 
descent methods or the Metropolis algorithms at zero temperature. The basis number TV for 
the periodic cell is varied from 1 to 16 to allow the sampling of variable-basis crystals over 
the entire pressure range defined by p m in and p max - If we find a lower enthalpy configuration 
than that for the diamond crystal, we add that configuration to the competitor list and 
repeat the two-step procedure. If we find no other states with lower enthalpy, we terminate 
the procedure. Using the final set of competitors, we have access to a large family of po- 
tentials that stabilize the diamond crystal with nearly-optimal pressure ranges. This allows 
us to choose a potential that has other useful qualities at the cost of only a small decrease 
in the pressure range over which the potential's ground state is the diamond crystal. One 
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such desirable property is for the lowest modes in the phonon spectrum to have the highest 
possible energies (i.e., the system is relatively stiff mechanically). 

After optimization under the aforementioned constraints (e.g., lack of floppy modes and 
convexity), we obtain the following optimized parameters for the potential function flTJ: 

at = -1.42324, a 2 = 0.713012. (3) 

Henceforth, we will refer to the potential function (CQ) with parameters defined by ([3]) as the 



diamond-1 or Dl potential;— see Fig. 1(a) The diamond crystal is the ground state of the Dl 
potential from p* = 0.0554 to p* = 0.1010. At these pressures, the corresponding densities 
and nearest-neighbor distances are, respectively, p* = 0.235 and p* = 0.303, and r^N = 
1.403a and r^N = 1.29a. The magnitudes of the second derivative at different distances 
have two simultaneous objects: (a) to stabilize the low-coordinated target structure, and 
(b) to discriminate against all competitors. Specifically, using the generalized coordination 
function formalism described in Ref. O, we determined that to obtain low-coordinated ground 
states with monotonic convex potentials, the magnitude of the second derivative must be 
large both below and near the nearest-neighbor distance of the low-coordinated ground state. 
It also must be small up to to the nearest-neighbor distance of the close-packed crystals and 
large up to the next-nearest-neighbor distance of the targeted low-coordinated ground state, 
after which it goes to zero. 




(a) 



(b) 



FIG. 1: (a) Optimized monotonic pair potential v(r) from Eq. (fTJ) using the parameters from ([3]): 
the Dl potential. |(b)| Second derivative d 2 v/dr 2 of the same optimized pair potential versus the 
distance r. 
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We also have attempted to use the same procedure to find a monotonic potential of the 
form ([1]) that has the closely related tetrahedrally-coordinated wurtzite crystal as its ground 
state. However, we have found that such a potential does not exist. This can be readily 
be explained by the fact that wurtzite and diamond crystal have the same short-range 
interparticle distances (first and second coordination shells), so the only way a potential of 
the form ([1]) can distinguish between the two is from its longer-range behavior. However, 
potentials of the form (JTJ decreases very rapidly at such distances, since they are dominated 
by the Gaussian factor. Therefore, the wurtzite crystal for the potential ([1]) has higher 
energy than the diamond crystal due to its third nearest neighbors being slightly closer. 

1 .2 i i i i i i 




Wavevector 



FIG. 2: Phonon spectrum in reduced units of the Dl potential for the diamond crystal at dimen- 
sionless pressure p* = 0.078 and density p* = 0.271. Only a representative subset of wavevectors 
that lie on paths connecting high-symmetry points (r, K, W, X, and Li) of the Brillouin zonei^ 
is shown. The Dl potential is chosen such that the ratio between the highest and lowest phonon 
frequencies at the X point is maximized. 

Figure [2] shows the phonon spectrum (which shows the strengths of the restoring forces for 
deformations associated with given wavevectors) of the Dl potential for the diamond crystal 
at a pressure in the middle of its stability range. Full phonon calculations (not only for the 
wavevectors shown in Fig. |2J) indicate that the diamond crystal is indeed mechanically stable 
for the Dl potential. The optimized Dl potential is selected among those potentials that 
produce a nearly optimal pressure range for the diamond ground state in order to achieve a 
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large ratio between the highest and lowest phonon frequencies at the X point. While only 
the X point was considered in the optimization, the Dl potential is also optimal for the 
other wavevectors that we have examined. 
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FIG. 3: Ground states of the Dl potential for a range of dimensionless pressures obtained from 
steepest descent for a basis up to N = 16. The crystal phases indicated from zero pressure to higher 
pressures are the 12-coordinated face-centered cubic (gray), the 8-coordinated body-centered cubic 
(cyan), a 2-coordinated hexagonal (orange), a 3-coordinated buckled rhombohedral graphite (blue), 
the 4-coordinated diamond (red), a 5/6-coordinated deformed diamond (green), a 6-coordinated 
buckled hexagonal (violet), and a 8-coordinated flattened-hexagonal closed-packed (yellow). Bonds 
are drawn between particles that are either nearest neighbors or nearly so as to guide the eye. 



The stable phases for the Dl potential at various pressures outside those for the diamond 
stability range are shown in Fig. [3j The phases are determined by repeatedly cooling disor- 
dered configurations at constant pressure using the aforementioned the variable-box energy 
minimization techniques and retaining the resultant configurations whose enthalpy are the 
lowest. We find that the diamond crystal is stable for 0.0554 < p* < 0.1010. Four neighbor- 
ing phases are particularly interesting. At low pressure, 0.0272 < p* < 0.0530, a hexagonal 
crystal phase, where the distance between hexagonal planes is shorter than the distance 
between particles from the same plane, is stable. This crystal has an effective coordination 
number of two, since the potential favors low-coordinated configurations over close-packed 
ones as long as the nearest-neighbor distance is not too small. Between this phase and the 
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diamond phase (0.0530 < p* < 0.0554), a low-coordinated rhombohedral graphite crystalline 
phase is stable. The rhombohedral graphite crystal is composed of stacked honeycomb crys- 
tals where each successive honeycomb layer is shifted in the same direction relative to the 
layer immediately below it (unlike standard graphite where the shift direction alternates 
between layers). The distance between the planes is 1.5 ~ 1.6 times larger than the nearest- 
neighbor distance within a layer, which is much less than that for actual graphite. At high 
pressure in the range 0.1155 < p* < 0.1315, the opposite happens, since the stable phase is 
a buckled simple hexagonal crystalline, where the nearest-neighbor distance within a hexag- 
onal plane is shorter than that between planes, resulting in a coordination number of six. 
Unlike the low-pressure hexagonal crystalline phase, this phase shows buckling: particles 
in the same layers are not perfectly aligned, but the distance between nearest neighbors 
stays constant. The transition between the high-pressure buckled hexagonal phase and the 
diamond phase (0.1010 < p* < 0.1155) consists of a highly deformed diamond crystal, for 
which particles have variable coordination numbers of either 5 or 6. The highest-pressure 
phase shown in Fig. [3] is a flattened- hexagonal closed-packed crystal, consisting of a distorted 
hexagonal closed-packed crystal with layer nearest-neighbor distances that are larger than 
the interlayer nearest-neighbor distance. This is not the stable phase for all p* > 0.1315; 
other phases arise at higher pressures. 



N 


Dl potential 


Star-polymer potential 


2 


96.89% 


91.41% 


4 


89.71% 


77.38% 


8 


62.13% 


54.28% 


16 


46.32% 


26.55% 


32 


24.57% 


8.93% 


64 


5.27% 


0.30% 



TABLE I: Frequency with which the ground-state diamond crystal is obtained from a steepest de- 
scent starting from a random configurations of iV particles. For each N, the frequency is calculated 
using 10000 trials, which results in standard deviations smaller than 0.5%. The Dl potential trials 
are carried out at p* = 0.078, while the star-polymer potential trials used p/(jgkBTf 3 ' 2 ) = 3.332 
(for which the ground state has a "packing fraction" r\ = 1.2) and an arm number / = 64.— 
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Configurational coordinate Configurational coordinate Configurational coordinate 

(a) (b) (c) 

FIG. 4: Schematic representations of three different types of energy landscapes as a function of the 
configurational coordinate. The boundaries of the basin of attraction associated with the global 



minima are indicated by the dashed vertical lines, (a) A relatively rough energy landscape. |(b)| An 



energy landscape with a deep and narrow global minimum, (c) An energy landscape with a broad 
and smooth global minimum. 

We employed the same rapid cooling method used to map out the phase diagram to 
quantify how easy it is for the system to reach the ground state. Table U compares the 
frequency with which the diamond crystal ground state is obtained using the Dl and the star- 
polymer— potentials for various numbers of particles, demonstrating an advantage of the Dl 
potential. The high frequency with which the Dl potential results in the diamond crystal is 
evidence that its energy landscape is smooth and possesses a broad global minimum, similar 



to the schematic illustration in Fig. 4(c) The decreasing frequencies for larger N are the 



consequence of the relative crudeness of our instant cooling method, which is ineffective at 
resolving large-scale defects. However, it is all the more remarkable that this method is 
capable of reaching the ground state with reasonable frequency despite using large bases, 
as opposed to, for example, a carefully-tuned simulated annealing procedure. Nevertheless, 
we have verified using simulated annealing techniques with a 256-particle system that the 
diamond crystal emerges as the ground state for the Dl potential. 

The work reported in this paper provides yet another example of the "inverse" statistical 
mechanical method to identify an appropriate interaction potential whose non-degenerate 
classical ground state is a pre-selected crystal structure. In general it is not guaranteed 
at the outset that such a targeted requirement has a solution. But in the present case of 
the fourfold coordinated diamond lattice previous studies have indeed established that this 
can be accomplished with pairwise additive isotropic potentials.- 1 ^ 1 ^ The existence of these 
examples establish that an infinite family of such interactions will produce the diamond 
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structure as its ground state, each member within some pressure (i.e., density) range. 

However merely stabilizing a given target structure is typically only part of the technical 
objective. In addition there tend to be properties that one wishes simultaneously to satisfy 
or to optimize. For the present study these have included maximizing the pressure range 
of ground-state stability for the diamond lattice, constraining the potential to monotonicity 
and convexity, maximizing the ratio of transverse to longitudinal acoustic sound speeds, and 
optimizing capture probability in the desired crystal basin from random initial configura- 
tions. It should be emphasized that the choice of such constraints and/or optimizations 
is not unique, but is driven by overall scientific objectives. Distinct choices obviously will 
identify distinct optimizing potential functions. 

The success at constructing diamond potentials naturally raises the question of whether 
the structure of that other macroscopic crystalline form of elemental carbon, graphite, might 
analogously be the classical ground state of an isotropic pair potential. At first sight this 
might seem easy, given the existence of potentials that generate the two dimensional analog, 
the honeycomb crystal.- However, the layered structure of this three-dimensional graphite 
allotrope, the stable form of elemental carbon at ambient conditions, with rather large 
interlayer separation and interlayer relative shift, realistically appears to present a formidable 
challenge. It may turn out that the graphite structure is unattainable with monotonically 
decreasing pair interactions. An ideal circumstance would be to find a potential whose 
classical ground state includes both the graphite structure (at low pressure) and the diamond 
structure (at elevated pressure), thus emulating reality. This ambitious joint requirement 
might require at least a combination of two-body and three-body interactions, suggesting a 
direction for future research. 

At present there is no known constraint on the complexity (basis of the unit cell) of a 
single-species target crystal structure that might be stabilized by an isotropic pair potential. 
But as the unit cell of a target crystal structure increases in size and geometric detail, it is 
reasonable to suppose that stabilizing isotropic pair potentials, if they exist, will necessarily 
also have to increase in range and complexity. Establishing such a connection constitutes 
another direction in which future studies should be focused. 
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